# ==============================================================================
# Title: strtools-0.1.R
# Description:
# Author:  Pierre Bady <pierre.bady@unil.ch>
# Date : 2010-09-10
# Version: 0.1
# Revision: Jan 26, 2011 9:05:22 AM
# Comments: RAS
# License: GPL version 2 or newer
# Copyright (C) 2011  Pierre Bady
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 2
# of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
# ==============================================================================

#--------------------------
# additional functions
# 2010-09-10 by pbady
#--------------------------
REbracketout <-function(x,replaceby="",...){
	x <- unlist(sapply(x,function(z) gsub("([(])|([)])",replaceby,z)))
	return(x)		
}
REzeroout <- function(x,replaceby="",...){
	x <- unlist(sapply(x,function(z) gsub("?([.][0])",replaceby,z)))
	return(x)	
}
REspaceout <- function(x,replaceby="",...){
	x <- unlist(sapply(x,function(z) gsub(" ",replaceby,z)))
	return(x)	
}
recode.STR2col <- function(df,bloc,sep="[-]",method="fuzzy",...){
# 2 methods , fuzzy or binary	
	if(!inherits(df,"data.frame"))
		stop("non convenient argument!")
	if(ncol(df)!=length(bloc))
		stop("non equal dimension")
	transfac <- function(x,level,type,sep=""){
		fx0 <- rep(0,length(level))
		fx1 <- factor(unlist(strsplit(x,sep)),level=level)
		if(type=="binary"){
			if(!all(is.na(fx1)))
				fx0[as.numeric(na.omit(fx1))] <- rep(1,length(na.omit(fx1)))
		}else{
			fx0 <- as.vector(table(fx1))
		}
		names(fx0) <- level
		fx0
	}
	bloc <- as.factor(bloc)
	df <- as.data.frame(apply(df,2,as.character),row.names=rownames(df),stringsAsFactors=FALSE)
	lw <- list()
	for(i in 1:nlevels(bloc)){
		auxi <- df[,bloc==levels(bloc)[i]]
		level1 <- na.omit(as.character(unlist(auxi)))
		level1 <- unlist(strsplit(level1,sep))
		level1 <- unique(level1)
		lw[[i]] <- apply(auxi,1,function(x) transfac(as.character(x),level1,method,sep))
	}
	names(lw) <- levels(bloc)
	w <- as.data.frame(t(do.call("rbind",lw)))
	nbloc <- unlist(lapply(lw,function(x) ifelse(is.null(nrow(x)),1,nrow(x))))
	auxinames <- NULL
	k2 <- 0
	for(i in 1:nlevels(bloc)){
		k1 <- k2 + 1
		k2 <- k2 + nbloc[i]
		if(nbloc[i]==1){
			tmp <- paste(levels(bloc)[i],unique(df[,bloc==levels(bloc)[i]][,1]),sep=".")
		}else{
			tmp <- paste(levels(bloc)[i],colnames(w)[k1:k2],sep=".")
		}
		auxinames <- c(auxinames,tmp)
	}
	colnames(w) <- auxinames
	attr(w,"bloc") <- nbloc
	attr(w,"levels") <- unlist(lapply(lw,rownames))
	attr(w,"fac") <- as.factor(rep(names(nbloc),nbloc))
	class(w) <- c("genocell",class(w))
	return(w)	
}
###### à verifier !!! pour les monoalléliques #####
recode.STR1col <- function(df,sep=",",method="fuzzy",...){
	if(!inherits(df,"data.frame"))
		stop("non convenient argument!")
	transfac <- function(x,level,type,sep=""){
		fx0 <- rep(0,length(level))
		fx1 <- factor(unlist(strsplit(x,sep)),level=level)
		if(type=="binary"){
			if(!all(is.na(fx1)))
				fx0[as.numeric(na.omit(fx1))] <- rep(1,length(na.omit(fx1)))
		}else{
			fx0 <- as.vector(table(fx1))
		}
		names(fx0) <- level
		fx0
	}
	df <- as.data.frame(apply(df,2,as.character),row.names=rownames(df),stringsAsFactors=FALSE)
	lw <- list()
	for(i in 1:ncol(df)){
		auxi <- df[,i]
		level1 <- na.omit(as.character(auxi))
		level1 <- unique(unlist(strsplit(level1,sep)))
		auxi <- t(sapply(auxi,function(x) transfac(as.character(x),level1,method,sep)))
		rownames(auxi) <- rownames(df) 
		lw[[i]] <- auxi
	}
	names(lw) <- colnames(df)
	w <- as.data.frame(do.call("cbind",lw))
	nbloc <- unlist(lapply(lw,function(x) ifelse(is.null(nrow(x)),1,ncol(x))))
	auxinames <- NULL
	k2 <- 0
	for(i in 1:length(nbloc)){
		k1 <- k2 + 1
		k2 <- k2 + nbloc[i]
		auxinames <- c(auxinames,paste(names(nbloc)[i],colnames(w)[k1:k2],sep="."))
	}
	colnames(w) <- auxinames
	attr(w,"bloc") <- nbloc
	attr(w,"levels") <- unlist(lapply(lw,colnames))
	attr(w,"fac") <- as.factor(rep(names(nbloc),nbloc))
	class(w) <- c("genocell",class(w))
	return(w)	
}
#### à modifier !!!! 2010-12-03 by pbady
merge.genocell <- function(x,y,...){
	if(!inherits(x,"genocell"))
		stop("non convenient argument")
	if(!inherits(y,"genocell"))
		stop("non convenient argument")
	vec <- intersect(names(attributes(x)$bloc),names(attributes(y)$bloc))
	lw <- list()
	xbloc <- attributes(x)$bloc
	xfac <- attributes(x)$fac
	ybloc <- attributes(y)$bloc
	yfac <- attributes(y)$fac
	ntot <- nrow(x)+nrow(y)
	for(i in 1:length(vec)){
		w1 <- x[,xfac==vec[i]]
		w2 <- y[,yfac==vec[i]]
		auxivec <- union(colnames(w1),colnames(w2))
		res <- matrix(0,nrow=ntot,ncol=length(auxivec))
		res[1:nrow(x),match(colnames(w1),auxivec)] <- as.matrix(w1)
		res[(nrow(x)+1):ntot,match(colnames(w2),auxivec)] <- as.matrix(w2)
		res <- as.data.frame(res)
		rownames(res) <- c(rownames(x),rownames(y))
		colnames(res) <- auxivec
		lw[[i]] <- res
		res <- NULL
	}
	names(lw) <- vec
	w <- as.data.frame(do.call("cbind",lw))
	nbloc <- unlist(lapply(lw,ncol))
	colnames(w) <- unlist(lapply(lw,colnames))
	attr(w,"bloc") <- nbloc
	attr(w,"levels") <- unlist(lapply(lw,colnames))
	attr(w,"fac") <- as.factor(rep(names(nbloc),nbloc))
	attr(w,"indica") <- c(rep(1,nrow(x)),rep(2,nrow(y)))
	class(w) <- c("genocell",class(w))
	return(w)	
}
prep.fuzzy.genocell <- function(x,...){
	if(!inherits(x,"genocell"))
		stop("non convenient argument")
	prep.fuzzy.var(x,attributes(x)$bloc,...)
}
heterozygosity <- function(x,...){
	UseMethod("heterozygosity")	
}
heterozygosity.genocell <- function(x,summary=TRUE,...){
	if(!inherits(x,"genocell"))
		stop("non convenient argument")
	fac <- attributes(x)$fac
	res <- list()
	for(i in 1:nlevels(fac)){
		w <-x[,fac==levels(fac)[i]]
		# ok si binaire, pas vrai si trisomy!
		homo <- unlist(apply(w,1,function(z) sum(z>0)==1 & !any(is.na(z))))
		Na <- unlist(apply(w,1,function(z) sum(z)==0 & !any(is.na(z))))
		hetero <- !homo & !Na
		mt2 <- unlist(apply(w,1,function(z) sum(z>0)>2 & !any(is.na(z))))
		res[[i]] <- cbind(homo,hetero,Na,mt2)
	}
	names(res) <- levels(fac)
	if(summary){
		res <- data.frame(do.call("rbind",lapply(res,function(z) apply(z,2, sum,na.rm=TRUE))))
		res$n <- apply(res[,1:2],1,sum,na.rm=TRUE)
		res$N <- apply(res[,1:3],1,sum,na.rm=TRUE)
		res$Nallele <- table(fac)
	}	
	return(res)
}
prep.DSMZ <- function(...){
	.NotYetImplemented()	
}
matchgenocell <- function(x,...){
	if(!inherits(x,"list)"))
		stop("non convenient argument!")
	.NotYetImplemented()
}
#----------------------------------------------------------------------------
# Evaluation values (see http://www.dsmz.de/fp/cgi-bin/help.html)
# based on dist.binary functions from R package ade4 (chessel et al. 2004
# if y is null, => dist.binary, method=5
# 2010-09-13 by pbady
#----------------------------------------------------------------------------
EVindex <- function(x,y=NULL,...){
	if (!(inherits(x, "data.frame") | inherits(x, "matrix"))) 
		stop("x is not a data.frame or a matrix")
	if(!is.null(y))
		if (!(inherits(y, "data.frame") | inherits(y, "matrix"))) 
			stop("y is not a data.frame or a matrix")
	x <- as.matrix(x)
	y <- as.matrix(y)	
	if (!is.numeric(x)) 
		stop("x must contain  numeric values")
	if (!is.numeric(y)) 
		stop("y must contain  numeric values")	
	if (any(x < 0)) 
		stop("non negative value expected in df")
	if (any(x < 0)) 
		stop("non negative value expected in df")	
	x <- as.matrix(1 * (x > 0))
	y <- as.matrix(1 * (y > 0))
	f1 <- function(df){
		a <- df %*% t(df)
		b <- df %*% (1 - t(df))
		c <- (1 - df) %*% t(df)
		d <- ncol(df) - a - b - c
		2 * a/(2 * a + b + c)
	}
	f2 <- function(df1,df2){
		a <- df1 %*% t(df2)
		b <- df1 %*% (1 - t(df2))
		c <- (1 - df1) %*% t(df2)
		d <- (1 - df1) %*% (1 - t(df2))
		2 * a/(2 * a + b + c)
	} 
	if(is.null(y)){
		w <- f1(x)
		w <- as.data.frame(w)
		if(!is.null(rownames(x))){
			rownames(w) <- rownames(x)
			colnames(w) <-  rownames(x)
		}
	}
	else{
		w <- f2(x,y)
		w <- as.data.frame(w)
		if(!is.null(rownames(x))) rownames(w) <- rownames(x)
		if(!is.null(rownames(y))) colnames(w) <- rownames(y)		
	}
	return(w)
} 
########## 2011-01-21 ########
simAIJ.marker <-function(marker=8,R=10){
	N <- marker*2
	a <- NULL
	for(i in 1:R){	
		x1 <- sample(marker:N,1)
		x2 <- sample(marker:N,1)
		a <- c(a,sum(sample(c(0,1),min(c(x2,x1)),replace=TRUE)))
	}
	return(a)
}
simAIJ.peak <-function(peak,R=10){
	a <- NULL
	for(i in 1:R){	
		a <- c(a,sum(sample(c(0,1),peak,replace=TRUE)))
	}
	return(a)
}
simSIM.marker <-function(marker=9,R=999){
	N <- marker*2
	w <- NULL
	for(i in 1:R){	
		x1 <- sample(marker:N,1)
		x2 <- sample(marker:N,1)
		a <- sum(sample(c(0,1),min(c(x2,x1)),replace=TRUE))
		w <- c(w,2*a/(x1+x2))
	}
	return(w)
}
simSIM2.marker <-function(marker=9,R=999){
	N <- marker*2
	w <- NULL
	for(i in 1:R){	
		x1 <- sample(marker:N,1)
		x2 <- sample(marker:N,1)
		a <- sum(sample(c(0,1),min(c(x2,x1)),replace=TRUE))
		w <- c(w,a/(x1+x2-a))
	}
	return(w)
}
simSIM3.marker <-function(marker=9,R=999){
	N <- marker*2
	w <- NULL
	for(i in 1:R){	
		x1 <- sample(marker:N,1)
		x2 <- sample(marker:N,1)
		a <- sum(sample(c(0,1),min(c(x2,x1)),replace=TRUE))
		w <- c(w,a)
	}
	return(w)
}
simSIM4.marker <- function(marker=8,R=9999){
	N <- 2*marker
	x1 <- sample(marker:N,1)
	x2 <- sample(marker:N,1)
	obs <- sum(sample(c(0,1),min(c(x2,x1)),replace=TRUE))
	sim <- simAIJ.peak(min(sum(x1),sum(x2)),R)
	w <- as.rtest(obs = obs, sim = sim, call = match.call())
	return(w)
}
simSIM4plus.marker <- function(marker=8,R=999,nrep=999){
	w <- NULL
	for(i in 1:R){
		w <- c(w,simSIM4.marker(marker=marker,R=nrep)$pvalue)
	}
	return(w)
}	
############### END! ###################	
